
#include <HYPRE.h>
#include <HYPRE_IJ_mv.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>

#include "matrices.h"
#include "errlog.h"
#include "glodef.h"

/****************************************************************************
 * Sparse_IJ  - compressed sparse column format
 ****************************************************************************/

typedef struct _SPARSE_IJ
{
    char    mtrtyp[3];   /* must be "ij" */
    int     neqns;
    HYPRE_IJMatrix A;
    
    int     *colstr;
    int     **aux_data;
    struct
    {
        unsigned int   structure : 1;
        unsigned int   assemble  : 1;
        unsigned int   changed   : 1;
    } flags;
} Sparse_IJ;

static void BeginAssembleSparse_IJ(Sparse_IJ *data)
{
    data->flags.structure = 0;
    data->flags.assemble = 1;
    data->flags.changed = 0;
    HYPRE_IJMatrixInitialize(data->A);
}

static void EndAssembleSparse_IJ(Sparse_IJ *data)
{
    if( data->flags.changed == 1)
    {
        warning_msg("%s: structure of matrix changed during reassembling.\n", __func__);
    }
    HYPRE_IJMatrixAssemble(data->A);
    data->flags.changed = 0;
    data->flags.structure = 0;
    data->flags.assemble = 0;
}

static void BeginStructureSparse_IJ(Sparse_IJ *data)
{
    data->flags.structure = 1;
    data->flags.assemble = 0;
    data->flags.changed = 0;
    
    if((data->colstr=malloc((1+data->neqns)*sizeof(int )))==NULL)
        goto alloc_error;
    if((data->aux_data=calloc(data->neqns, sizeof(int  *)))==NULL)
        goto alloc_error;
    
    for(int  i=0; i<data->neqns; i++)
    {
        data->colstr[i] = 1;
        if((data->aux_data[i]=calloc(1, sizeof(int )))==NULL)
            goto alloc_error;
        data->aux_data[i][0]=i;
    }
    data->colstr[data->neqns] = 0;
    
    return;
    alloc_error:
        error_msg("%s: memory allocation error.\n", __func__);
}

static void EndStructureSparse_IJ(Sparse_IJ *data)
{
    int  i;
    data->flags.structure = 0;
    data->flags.assemble = 0;
    data->flags.changed = 0;
    
    HYPRE_IJMatrixSetRowSizes(data->A, data->colstr);

    free(data->colstr);
    data->colstr = NULL;
    
    for(i=0;i<data->neqns;i++)
        free(data->aux_data[i]);
    free(data->aux_data);
    data->aux_data = NULL;

    return;
}




static Sparse_IJ *AllocSparse_IJ(int  neqns, int sym, char *mtrtyp)
{
    /* here sym refers to matrix values.  If matrix is structurally symmetric,*/
    /* but unsymmetric values, then sym must be FALSE */
    Sparse_IJ *data;
    data = calloc(1, sizeof(Sparse_IJ));
    if(data==NULL)
        error_msg("Cannot allocate Sparse_Matrix in IJ storage.\n");
    data->neqns = neqns;
    data->mtrtyp[0]=tolower(mtrtyp[0]);
    data->mtrtyp[1]=tolower(mtrtyp[1]);
    data->mtrtyp[2]=0;
    if( strcmp(data->mtrtyp, "ij")!=0 )
        error_msg("Invalid IJ matrix type %s\n", data->mtrtyp);
    HYPRE_IJMatrixCreate(MPI_COMM_WORLD, 0, neqns-1, 0, neqns-1, &data->A);
    HYPRE_IJMatrixSetObjectType(data->A, HYPRE_PARCSR);
    return data;
}

static void CleanSparse_IJ(Sparse_IJ *data)
{
    HYPRE_IJMatrixInitialize(data->A);
}

static void DestroySparse_IJ(Sparse_IJ *data)
{
    if(data)
    {
        if(data->colstr) free(data->colstr);
        if(data->aux_data)
        {
            for(int i=0;i<data->neqns; i++)
                free(data->aux_data[i]);
            free(data->aux_data);
        }
        HYPRE_IJMatrixDestroy(data->A);
        free(data);
    }
}

static void PrintSparse_IJ(FILE *F, Sparse_IJ *data)
{
    error_msg("%s: operation is not available.\n", __func__);
}

/*
 * static void _SM_IncreaseRowArrays(Sparse_Row *R, int incr)
 * {
 * void *ptr;
 * int ns;
 *
 * ns = R->EntriesCount + incr;
 * ptr = realloc(R->Columns, ns*sizeof(int));
 * if( ptr == NULL )
 * error_msg("cannot allocate memory for sparse matrix.");
 * R->Columns=(int*)ptr;
 * ptr = realloc(R->Values, ns*sizeof(double));
 * if( ptr == NULL )
 * error_msg("cannot allocate memory for sparse matrix.");
 * R->Values=(double*)ptr;
 * R->EntriesCount = ns;
 * }
 */

static double* _SM_FindEntrySparse_IJ(Sparse_IJ *data,
        int  row, int  col, int add)
{
    int  i;
    void *ptr;
    
    if(data->mtrtyp[0]=='s' && (data->mtrtyp[1]=='p' || data->mtrtyp[1]=='s'))
    {
        /* storage is symmetric and values are symmetric */
        /* we only store the lower triangular part, so ...  */
        if(row < col)
        { /* swap row <-> col */ i = row; row = col; col = i; }
    }
    
    if(data->flags.structure)
    {
        for(i=0; i<data->colstr[col]; i++)
            if(row == (data->aux_data[col][i]))
                return data->values;
        int  foo = ++data->colstr[col];
        ptr = realloc(data->aux_data[col], foo*sizeof(int ));
        if(ptr==NULL) goto alloc_error;
        data->aux_data[col] = ptr;
        #ifdef  IJ_SORT_COLS
                for(i=foo-1;i>0;i--)
                {
            if(row < (data->aux_data[col][i-1]))
                data->aux_data[col][i] = data->aux_data[col][i-1];
            else
            {
                data->aux_data[col][i] = (row);
                break;
            }
                }
        if(i==0) data->aux_data[col][0] = (row);
        #else
                data->aux_data[col][foo-1] = (row);
        #endif
                return data->values;
    }
    
    /* TODO: wtf */
    
    for(i=f2c(data->colstr[col]); i<f2c(data->colstr[col+1]); i++)
    {
        if(row == f2c(data->rowind[i]))
            return data->values + i;
        #ifdef IJ_SORT_COLS
                if(row < f2c(data->rowind[i]))
                    break;
        #endif
    }
    
    if(!add) return NULL;
    
    data->flags.changed = 1;
    
    /* add entry (row,col) at position i in arrays, push next entries back */
    for(col++; col<=data->neqns; col++) data->colstr[col]++;
    
    int  nentries = f2c(data->colstr[data->neqns]);
    
    ptr=realloc(data->values, nentries*sizeof(double));
    if(ptr==NULL) goto alloc_error;
    data->values=ptr;
    
    ptr=realloc(data->rowind, nentries*sizeof(int ));
    if(ptr==NULL) goto alloc_error;
    data->rowind=ptr;
    
    memmove(data->values+i+1, data->values+i, (nentries-i-1)*sizeof(double));
    memmove(data->rowind+i+1, data->rowind+i, (nentries-i-1)*sizeof(int ));
    
    data->rowind[i]=c2f(row);
    data->values[i]=0.0;
    return data->values+i;
    
    alloc_error:
        error_msg("%s: cannot allocate memory for sparse matrix.", __func__);
}

static void _SM_MultV_Sparse_IJ(Sparse_IJ *data,
        const double *vect, double **res)
{
    int  nentries = f2c(data->colstr[data->neqns]);
    int  i, row, col;
    double val;
    col=0;
    memset(*res, 0, data->neqns*sizeof(double));
    if(((strcmp(data->mtrtyp, "ss")==0) || (strcmp(data->mtrtyp, "sp")==0)))
        for(i=0;i<nentries;i++)
        {
        while(i>=f2c(data->colstr[col+1])) col++;
        row=f2c(data->rowind[i]);
        val=data->values[i];
        (*res)[row] += val * vect[col];
        if(row!=col) (*res)[col] += val * vect[row];
        }
    else
        for(i=0;i<nentries;i++)
        {
        while(i>=f2c(data->colstr[col+1])) col++;
        row=f2c(data->rowind[i]);
        val=data->values[i];
        (*res)[row] += val * vect[col];
        }
}

static void _SM_MultTransV_Sparse_IJ(Sparse_IJ *data,
        const double *vect, double **res)
{
    int  nentries = f2c(data->colstr[data->neqns]);
    int  i, row, col;
    double val;
    col=0;
    memset(*res, 0, data->neqns*sizeof(double));
    if(((strcmp(data->mtrtyp, "ss")==0) || (strcmp(data->mtrtyp, "sp")==0)))
        for(i=0;i<nentries;i++)     /* symmetric is the same as straight */
        {
        while(i>=f2c(data->colstr[col+1])) col++;
        row=f2c(data->rowind[i]);
        val=data->values[i];
        (*res)[row] += val * vect[col];
        if(row!=col) (*res)[col] += val * vect[row];
        }
    else
        for(i=0;i<nentries;i++)     /* transpose multiplication */
        {
        while(i>=f2c(data->colstr[col+1])) col++;
        row=f2c(data->rowind[i]);
        val=data->values[i];
        (*res)[col] += val * vect[row];
        }
}

/*****************************************************************************/
